!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_initial (base)
!==================================================================================================================================|

   use parameters
   use block_define

   implicit none

   integer(4) :: i, i_block, i_total, i_neigh, i_level
   integer(4) :: zerox

   real(8) :: vx, vy, vz, b2, v2, c2, s2, ca2, cfx, pr

   type(block), pointer, intent(inout) :: base

   type(block), pointer :: child

!----------------------------------------------------------------------------------------------------------------------------------|
   zerox = nx_total * zero_posx / rangex + 1
   base%dx = rangex / dble(nx_total - 2 * ng)

   base%x(zerox) = base%dx * 0.5d0
   
   do i = zerox + 1, nx_total
      base%x(i) = base%x(i-1) + base%dx
   enddo
   do i = zerox-1, 1, -1
      base%x(i) = base%x(i+1) - base%dx
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!   do i = 1, nx_total
!      pr = 1.0d0
!      base%ro(i) = 0.25d0
!      if (base%x(i) .ge. -0.40 .and. base%x(i) .le. -0.20d0) then
!         base%ro(i) = 0.5d0
!      endif
!      if (base%x(i) .ge. 0.0 .and. base%x(i) .le. 0.2d0) then
!         base%ro(i) = 0.25 * sin(5.0d0 * pr * base%x(i)) + 0.25d0
!      endif
!      base%mx(i) = 1.0d0 * base%ro(i)
!      base%my(i) = 0.0d0
!      base%mz(i) = 0.0d0
!      base%bx(i) = 0.0d0
!      base%by(i) = 0.0d0
!      base%bz(i) = 0.0d0
!      base%en(i) = pr / gmm1 + base%mx(i) * base%mx(i) / base%ro(i) / 2.0d0
!   enddo

!   do i = 1, nx_total
!      pr = 1.0d0
!      base%ro(i) = (sin(2.0d0 * pi * base%x(i)) + 2.0d0) / 3.0d0
!      base%mx(i) = 1.0d0 * base%ro(i)
!      base%my(i) = 0.0d0
!      base%mz(i) = 0.0d0
!      base%bx(i) = 0.0d0
!      base%by(i) = 0.0d0
!      base%bz(i) = 0.0d0
!      base%en(i) = pr / gmm1 + base%mx(i) * base%mx(i) / base%ro(i) / 2.0d0
!   enddo

!   do i = 1, nx_total
!      if (base%x(i) .lt. 0) then
!         base%ro(i) = 0.445d0
!         base%mx(i) = 0.311d0
!         base%my(i) = 0.0d0
!         base%mz(i) = 0.0d0
!         base%bx(i) = 0.0d0
!         base%by(i) = 0.0d0
!         base%bz(i) = 0.0d0
!         base%en(i) = 8.928d0
!      else
!         base%ro(i) = 0.500d0
!         base%mx(i) = 0.000d0
!         base%my(i) = 0.0d0
!         base%mz(i) = 0.0d0
!         base%bx(i) = 0.0d0
!         base%by(i) = 0.0d0
!         base%bz(i) = 0.0d0
!         base%en(i) = 1.4275d0
!      endif
!   enddo

!   do i = 1, nx_total
!      if (base%x(i) .lt. 0) then
!         base%ro(i) = 1.0d0
!         base%mx(i) = 0.0d0
!         base%my(i) = 0.0d0
!         base%mz(i) = 0.0d0
!         base%bx(i) = 0.0d0
!         base%by(i) = 0.0d0
!         base%bz(i) = 0.0d0
!         pr = 1.0d0
!      else
!         base%ro(i) = 0.125d0
!         base%mx(i) = 0.0d0
!         base%my(i) = 0.0d0
!         base%mz(i) = 0.0d0
!         base%bx(i) = 0.0d0
!         base%by(i) = 0.0d0
!         base%bz(i) = 0.0d0
!         pr = 0.1d0
!      endif
!      base%en(i) = pr/gmm1 + base%mx(i) * base%mx(i) / base%ro(i) / 2.0d0
!   enddo

!   do i = 1, nx_total
!      if (base%x(i) .lt. 0) then
!         base%ro(i) = 1.0d0
!         base%mx(i) = 0.0d0
!         base%my(i) = 0.0d0
!         base%mz(i) = 0.0d0
!         base%bx(i) = 0.0d0
!         base%by(i) = 0.0d0
!         base%bz(i) = 0.0d0
!      else
!         base%ro(i) = 0.125d0
!         base%mx(i) = 0.0d0
!         base%my(i) = 0.0d0
!         base%mz(i) = 0.0d0
!         base%bx(i) = 0.0d0
!         base%by(i) = 0.0d0
!         base%bz(i) = 0.0d0
!      endif
!      pr = 1.0d0
!      base%en(i) = pr/gmm1 + base%mx(i) * base%mx(i) / base%ro(i) / 2.0d0
!   enddo

   do i = 1, nx_total
      if (base%x(i) .lt. 0) then
         base%ro(i) = 1.0d0
         base%mx(i) = 0.0d0
         base%my(i) = 0.0d0
         base%mz(i) = 0.0d0
         base%bx(i) = 0.75d0
         base%by(i) = 1.0d0
         base%bz(i) = 0.0d0
         pr = 1.0d0
      else
         base%ro(i) = 0.125d0
         base%mx(i) = 0.0d0
         base%my(i) = 0.0d0
         base%mz(i) = 0.0d0
         base%bx(i) = 0.75d0
         base%by(i) = -1.0d0
         base%bz(i) = 0.0d0
         pr = 0.1d0
      endif
      base%en(i) = pr/gmm1 + base%mx(i) * base%mx(i) / base%ro(i) / 2.0d0 +                                                        &
         (base%bx(i) * base%bx(i) + base%by(i) * base%by(i) + base%bz(i) * base%bz(i)) / 2.0d0
   enddo

   max_speed = 0.d0
   do i = 1, nx_total
      vx = base%mx(i) / base%ro(i)
      vy = base%my(i) / base%ro(i)
      vz = base%mz(i) / base%ro(i)
      b2 = base%bx(i) * base%bx(i) + base%by(i) * base%by(i) + base%bz(i) * base%bz(i)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (base%en(i) - v2 * base%ro(i) / 2.0d0 - b2 / 2.0d0) * gmm1
      c2 = gm * pr / base%ro(i)
      s2 = c2 + b2 / base%ro(i)
      ca2 = base%bx(i) * base%bx(i) / base%ro(i)
      cfx = sqrt ((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / 2.0d0)
      max_speed = max (max_speed, cfx + abs(vx))
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!     initial blocks
!----------------------------------------------------------------------------------------------------------------------------------|
   write (base%id, "('B', i9.9)") id_index
   base%lv = 0
   base%nc = ini_blocks
   base%nx = nx_total
   base_dt = cfl * base%dx / max_speed
   base%dt = base_dt
   base%t = 0.0d0

   id_index = id_index + 1

   do i_block = 1, base%nc
      child => base%child(i_block)

      child%lv = 1
      child%nx = nx_block
      child%p_pos = i_block
      child%dx = base%dx
      child%dt = base_dt
      child%t = base%t
      child%pp => base

      allocate (child%x(child%nx), child%ro(child%nx), child%mx(child%nx), child%my(child%nx), child%mz(child%nx),                 &
         child%bx(child%nx), child%by(child%nx), child%bz(child%nx), child%en(child%nx),                                           &
         child%flux_bnd(n_neighs, sys), child%neigh(n_neighs, n_levels))

      write (child%id, "('B', i9.9)") id_index

      id_index = id_index + 1

      do i = 1, child%nx
         i_total = (i_block - 1) * nx_no_bnd + i
         child%x(i) = base%x(i_total)
         child%ro(i) = base%ro(i_total)
         child%mx(i) = base%mx(i_total)
         child%my(i) = base%my(i_total)
         child%mz(i) = base%mz(i_total)
         child%bx(i) = base%bx(i_total)
         child%by(i) = base%by(i_total)
         child%bz(i) = base%bz(i_total)
         child%en(i) = base%en(i_total)
         if (isnan(child%bx(i))) then
            print*, child%bx(i), i
         endif
      enddo

      do i_neigh = 1, n_neighs
      do i_level = 1, n_levels
         child%neigh(i_neigh, i_level)%pblock => null()
      enddo
      enddo

!     free boundary
      if (i_block .eq. 1) then
         if (ini_blocks .gt. 1) child%neigh_lv(2) = 1
         if (ini_blocks .gt. 1) child%neigh(2, 1)%pblock => base%child(i_block + 1)
      elseif (i_block .eq. ini_blocks) then
         if (ini_blocks .gt. 1) child%neigh_lv(1) = 1
         if (ini_blocks .gt. 1) child%neigh(1, 1)%pblock => base%child(i_block - 1)
      else
         child%neigh_lv(1) = 1
         child%neigh_lv(2) = 1
         child%neigh(2, 1)%pblock => base%child(i_block + 1)
         child%neigh(1, 1)%pblock => base%child(i_block - 1)
      endif
 
!     periodic boundary
!      if (i_block .eq. 1) then
!         child%neigh_lv(1) = 1
!         child%neigh_lv(2) = 1
!         child%neigh(2, 1)%pblock => base%child(i_block + 1)
!         child%neigh(1, 1)%pblock => base%child(ini_blocks)
!      elseif (i_block .eq. ini_blocks) then
!         child%neigh_lv(1) = 1
!         child%neigh_lv(2) = 1
!         child%neigh(2, 1)%pblock => base%child(1)
!         child%neigh(1, 1)%pblock => base%child(i_block - 1)
!      else
!         child%neigh_lv(1) = 1
!         child%neigh_lv(2) = 1
!         child%neigh(2, 1)%pblock => base%child(i_block + 1)
!         child%neigh(1, 1)%pblock => base%child(i_block - 1)
!      endif
   enddo

   return
end subroutine model_initial
